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ACCURATE  EFFICIENT  EVALUATION  OF  CUMULATIVE  OR 
EXCEEDANCE  PROBABILITY  DISTRIBUTIONS  DIRECTLY  FROM 
CHARACTERISTIC  FUNCTIONS 
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INTRODUCTION 

The  performance  of  a  signal  processor  can  often  be  evaluated  in  terms  of 
the  characteristic  function  of  the  decision  variable,  either  numerically  or  in 
closed  form;  see  for  example,  refs.  1  and  2.  However,  a  closed  form  for  the 
corresponding  probability  density  function  or  cumulative  distribution  function 
is  seldom  available,  and  numerical  procedures  must  be  employed.  Several  such 
procedures  have  been  published  in  the  literature,  refs.  3-8.  However  they 
have  limited  accuracy  or  they  require  extensive  storage  or  analytical 
manipulations  and  calculations. 

We  present  a  technique  which  is  limited  in  accuracy  only  by  the  round-off 
noise  of  the  computer  or  by  the  errors  of  the  special  functions  requi-ed  in 
the  characteristic  function  calculation.  The  amount  of  storage  depends  only 
on  the  number  of  cumulative  or  exceedance  distribution  function  values 
requested  and  does  not  influence  the  accuracy  of  the  final  probability 
values.  The  entire  cumulative  and  exceedance  distribution  function  values 
result  as  the  output  of  one  fast  Fourier  transform  (FFT).  The  si2e  of  the  FFT 
dictates  the  storage  required  and  the  spacing  of  the  calculated  probability 
values,  but  not  their  accuracy. 

The  addition  and  subtraction  of  integrand  functions  given  in  ref.  7  can 
be  entirely  circumvented  and  yet  enable  use  of  an  FFT,  through  proper 
manipulation  of  the  origin  contribution  of  the  characteristic  function. 
Specific  connections  with  past  results  will  be  noted  at  appropriate  points  in 
the  derivations. 
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DERIVATION  OF  PROCEDURE 


Shifted  Random  Variable 


The  primary  random  variable  of  interest  is  the  real  quantity  x  with  given 
characteristic  function  fx(?)  which  is  related  to  the  probability  density 
function  px  of  random  variable  x  via  Fourier  transform  * 

fx(5 )  «Jdv  exp(i$v)  Px(v)*  (D 

We  define  secondary  random  variable  y  as 

y  a  x+b,  (2) 

where  bias  (shift)  b  is  a  constant,  chosen  such  that  random  variable  y  has 
insignificant  probability  of  being  less  than  zero.  However,  we  also  pick  b  as 
small  as  possible,  so  that  the  characteristic  function  of  y, 

f  (?)  »  f  (?)  exp(ibf),  (3) 

y  * 

will  vary  slowly  with  In  fact,  b  can  be  negative,  as  for  example  if  x  were 
limited  to  values  larger  than  some  positive  threshold.  The  approach  here  is 
not  limited  to  positive  random  variables  x,  as  were  some  of  the  results  in 
ref.  7,  but  is  applicable  to  any  random  variable  distribution. 

By  way  of  example,  for  an  exponential  probability  density  function  for 
random  variable  x,  we  choose  b*0;  while  for  a  Gaussian  random  variable, 
b-u  *8ax  yields  a  probability  less  than  IE-15  of  y  being  negative.  The 
probability  density  function  of  random  variable  y  therefore  appears  as 
depicted  in  figure  1. 


•  Integrals  and  sums  without  limits  are  over  (-*>,♦*>). 
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figure  1.  Probability  Density  Function  of  Secondary  Random  Variable  y 


The  cumulative  distribution  functions  of  random  variables  y  and  x  are 
related  according  to 
v 

Jdt  p  (t)  -  P  (v)  *  Px(v-b);  Pxfw)  =  Py(v-b).  (4) 

* —  G# 

Tnus  we  can  inspect  P  (v)  in  the  neighborhood  of  v»-b  (the  lower  edge  of 
interest  of  x)  by  looking  at  cumulative  distribution  function  Py(v)  in  the 
neighborhood  of  v=»0.  More  precisely,  we  will  investigate  P  (v)  for  values 
of  v  greater  than  zero,  since  this  is  the  region  of  significant  variation  of 
P  (v);  this  is  called  the  positive  neighborhood  of  v»0. 

Approximation  to  Cumulative  Distribution  Function 


From  ref.  4,  eq.  7,  we  have  the  cumulative  distribution  function  of 
random  variable  y  ir,  terms  of  the  characteristic  function  according  to 


-fa?  (5) 

O 

where  we  have  defined  auxiliary  function 


g(f.v)  o 


Id 


|exp(-‘fv) 


Observe  for  later  use  that 


U) 


[  t 

g(0*,v)  =•  l  im  Ia<(l-ifv)  - 

t**  . 


u  -v 


where  yy  is  the  mean  of  random  variable  y. 


(?) 

3 


-riftry- 


TR  No.  7023 


For  v  in  the  neighborhood  of  zero,  exp(-i$v)  in  (6)  varies  slowly  with  f, 
and  we  have  the  approximation,  via  the  Trapezoidal  rule,  to  (5)  as 

+«© 

Py(v)  *  J  -  f  9(0*,v)  -  a  g(na,v)  s  C(v),  (3) 

n=l 

where  the  right-hand  side  of  (8)  has  been  defined  as  C(v).  Here,  a  is  the 
sampling  interval  in  $,  and  is  small  enough  to  track  changes  in  exp(-ifv)* 
fy ($)/£.  We  choose  the  Trapezoidal  rule  in  (8)  over  other  integration  rules, 
such  as  Simpson's  rule,  because  it  results  in  minimum  aliasing  for  Fourier 
transforms  relative  to  all  other  rules;  see  appendix  A  for  elaboration  and 
proof . 

Observe  from  (8)  that 

Py(0)  *  C(0)  means  C{0)  a  0,  (9) 

since  ( 0 )  is  insignificant  by  the  choice  of  b  in  (2);  this  relation  will 

be  used  later. 


TR  No.  7023 


The  removal  of  the  imaginary  operation  from  within  the  summation  in  (10)  i ;  a 
crucial  step;  it  does  not  create  a  problem  in  divergence  since  n  >  0.  This  is 
in  contrast  with  the  integral  of  (5)  and  (6),  where  removal  of  the  imaginary 
operation  would  create  a  divergent  integral.  This  postponement  of  the  removal 
of  the  imaginary  ooeration,  until  after  the  approximation  to  the  integral  was 
developed  in  (8),  is  the  major  difference  with  the  results  in  ref.  7. 

Taking  a  derivative  of  (11),  we  obtain 

C  (V)  “  17  *  27  ^eXP^if,AV)  f  y  ( n  a )  = 

n^O 

a  ^  ^exp(-inav)  fy(na)  = 
n 

«  ^  exp(-ifv)  fy(J)  *£(?)• 

«  Py(v)«  S2s(v}  =  ^Py(V~n  f1)  S  PylV).  (12) 

n 

where  infinite  impulse  train 

-  ^l(y-na),  (13) 

n 

where®  denotes  convolution,  and  where  we  have  used  the  relation 

J*  dt  exp(-i«t)  a&A(t)  *^,t(w).  (W! 


This  last  result  follows  from  ref.  9,  p.  28,  rule  11,  with  yft)  -  S(t),  T„a, 
f«l/T,  and  W*f.  Relation  (12)  indicates  that  C’(v)  is  an  infinitely  aliased 
version  of  the  probability  density  function  py(v),  with  resultant  period 
t»/a  in  v.  For  small  enough  stapling  increment  a  in  (8),  there  will  be  very 
lutle  overlap  of  the  displaced  versions  of  p  in  (12),  thereby  yielding  the 
good  approximation 
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Relation  (16)  is  an  exact  relation,  snowing  that  C(v)  is  the  integral  of  the 
infinitely  aliased  version  of  p  (v),  starting  at  v-0,  plus  an  additive 
constant  which  is  substantially  zero;  see  (9). 


So  for  v  in  the  positive  neighborhood  of  zero,  (4),  (8),  (16),  and  (9) 

yield 


Px(v-b)  -  Py(v)  *  C( v )  -  C( 0)  *  Py(v)  a  Py(v).  (13) 

Thus  the  quantity  we  want,  the  left-most  term  in  (18),  is 
wei  1- approximated  by  calculated  quantity  C(v),  which  itself  is  approximately 
the  integral  of  the  infinitely  aliased  version  of  p  (v). 
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Calculation  of  C ( v ) 


2nk  . 


Let  v  =  -j^-  in  (10),  where  M  and  k  are  aroitrary  integers.  Then 


C(^r)  =  f  +  Im[2  exp(-i21 


nk/M)- 


fu(na) 


irh 


,n=l 


>  - 


f**  "J 

=  \  +  I  ~  «  Im<  exP(~i2ltnk/M)  znr> 

L  riaO  J 


wnere  we  define  complex  sequence 


(19) 


^  i  v  for  n=0 
|f^(na)/n  for  n>l 


Now  define  collapsed  sequence  (ref.  7,  pp.  13-16)  as 


+*> 


J»0 


Zn+Mj  for  °<n£ti-1- 


(20) 


(21) 


Tnen  since  zn  receives  tne  same  weight  as  zn+^j  <n  (19),  regardless  of  the 
value  of  k,  (19)  can  oe  expressed  as 


c's?i 


exp(-i2«nk/M) 


(22) 


Relation  (2z)  is  exact  and  valid  for  all  x.  Since  we  are  only  interested 
in  tne  positive  netgnoornood  of  v«0  in  (18;,  we  confine  attention  in  (22)  to 
U  <_  k  ^  M-i.*  Relation  1221  can  then  oe  accomplished  Dy  an  H-point  FFT  if  N 
is  chosen  to  De  a  power  of  2.  Notice  that  storage  only  for  the  N  complex 
numoers  lf5  (21)  is  required,  even  though  tne  {zn}  sequence  in  (20)  is 

of  infinite  length. 


*  Values  for  otner  k  are  avai ladle  from  (22)  when  we  ooserve  that 

‘Hr*1)  ■ 1  • c  (sf) for  a"  *• 
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Observe  that  the  size  of  M  in  no  way  affects  the  error  of  the  calculation  of 

2irk 

^Ma-^  or  estimati°n  of  Py(v).  Rather,  M  specifies  the  spacing  at  which  C(^~)  is 
calculated,  and  can  be  coarse  if  desired.  The  accuracy  of  the  estimate  of 
Py(v)  ^  governed  thus  far  by  a,  through  the  aliasing  depicted  in  figure  2. 

Reference  to  (18)  now  yields 

Px(ffiT  “  3  for  °<k<M-U  (2a) 

where  the  latter  quantity  is  given  by  (22).  Thus  the  M-point  FFT  sweeps  out 
the  argument  range  (-b,-b+2t/A)  for  the  cumulative  distribution  function  ?x. 

If  we  want  the  exceedance  distribution  function  of  y  instead  of  the 
cumulative  distribution  function,  we  use  (13)  and  (22)  to  get 


1  r / 2wk  v  1  k  ^  1 

1  '  “  1 


/M-l  -v 

“  j?  +  7  In>  exp(-i2*nk/M)  zi  for  0  <  k  <  M-l.  (24) 
n=0  J 


(By  the  footnote  to  (22),  we  have  1-C(2«/a)  -  -C(0).) 


Since  Uy  must  oe  Known  in  (20)  in  order  to  use  this  approach,  we  need 
tne  mean  Uj,  of  random  variable  x,  since  from  (2) 


U 


V 


ux4b. 


(25) 


The  quantity  u<  can  be  found  analytically  from  characteristic  function 
fx(J)  according  to 


f'(0)  »  i 
x 


(25) 


see  (1). 
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In  addition  to  the  error  caused  by  aliasing  associated  with  nonzero 
sampling  increment  a,  an  additional  error  occurs  because  we  cannot  calculate 
all  the  coefficients  fz^  in  (20)  and  (21)  out  to  n=+*.  Rather,  we 
terminate  the  calculation  at  integer  n=N,  such  that  JzJ  is  sufficiently 
small  as  to  be  negligible  for  n  _>  N.  Letting 

L  =  Na,  (27) 

this  is  equivalent  to  ignoring  the  contribution  to  (5)  of  the  tail  error 

+  00  +00 

r  r  fva) 

-  \d?  g(?,v)  =  -Im  \  d£  exp(-ijv)  .  (28) 

«  j  if  y 

L  L 

If  the  asymptotic  behavior  of  fy(£)  for  large  J  is  known,  this  error  can 
sometimes  be  evaluated  in  closed  form  and  used  to  ascertain  an  adequate  value 
rtf  L.  Instead,  we  hive  observed  that  tail  error  (28)  causes  a  characteristic 
low-level  sinusoidal  variatic-,1  in  the  calculated  cumulative  distribution 
function  for  small  v  near  0,  and  in  the  calculated  exceedance  distribution 
function  for  large  v  near  2i,/a.  When  this  sinusoidal  variation  is  deemed 
excessive,  L  can  be  increased  until  the  effect  disappears  or  decreases  to 
acceptable  levels.  This  tria1  and  error  approach  avoids  the  necessity  of 
analytically  upoer  -bounding  the  magnitude  of  error  (28),  which  is  often  very 
tedious  and  generally  nessimistic. 

So  there  are  two  errors  to  be  concerned  with:  aliasing  due  to  nonzero 
sampling  interval  A  and  tail  error  due  to  non  infinite  limit  L.  Later 
examples  will  demonstrate  how  t'.  ,-e  errors  manirest  themselves  in  the 
cumulative  and  exceedance  distribution  functions  and  how  they  can  be 
controlled  by  a  trial  and  error  approach. 
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Relation  to  Requicha's  Method,  ref.  5 


From  ref.  5,  eqs.  7,  9,  10,  the  cumulative  distribution  function  is  given 
by  an  expression  that  can  be  manipulated  into  the  form  (using  current  notation) 


fi/2 


Fk’R-7I“'f 


2  exp(-i2.kn/M)  ^  .  1  , 
n=l  (n=l 


(29) 


Although  this  is  similar  to  the  upper  line  of  (19)  here,  it  differs  in  several 
important  respects: 


1.  F^  does  not  use  mean  at  all;  it  is  therefore  not  using  a 
direct  approximation  to  the  specified  integral  in  (5)  and  (6). 

2.  From  (29),  there  follows  Fg  =  0,  FM  =  1;  however,  these  results 

are  not  strictly  true  for  the  actual  cumulative  distribution  function  at  these 
end  points,  thereby  leading  to  poor  estimates  in  the  neighborhoods  of  these 
points.  This  is  due  to  the  arbitrary  origin  established  in  ref.  5,  eq.  6. 

3.  The  sums  in  (29)  utilize  characteristic  function  samples  f  (na) 
only  for  n  <  M/2,  where  M  is  the  size  of  the  FFT.  This  is  a  very  severe  and 
unnecessary  restriction;  in  fact,  the  sum  on  n  in  (29)  ought  to  be  conducted 
to  the  point  where  the  tail  contribution,  (28),  is  negligible,  regardless  of 
the  value  of  M. 


4.  In  ref.  5,  if  eq.  4  is  substituted  into  eq.  1,  and  the  summation 
limits  are  extended  to±»,  we  get  exactly  the  second  line  of  (12)  here.  When 
the  probability  density  function  is  integrated  to  get  the  cumulative 
distribution  function  in  ref.  5,  eq.  6,  the  resultant  cumulative  distribution 
function  is  arbitrarily  set  to  zero  at  v=0.  We  instead  have  from  (9)  and  (17), 


which  is  small,  but  not  necessarily  zero.  This  consideration  is  very 

important  on  the  tails  of  the  cumulative  and  exceedance  distribution  functions. 
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EXAMPLES 

Programs  for  the  following  five  examples  are  listed  in  appendix  B. 
1.  Chi-Square 


A  chi-square  variate  of  2K  degrees  of  freedom  has  probability  density 
function  (ref.  10) 


p  (v)  -  vK-^-  e-xPl~v/2i-  for  v  >  0  (35) 

x  2K(K-l)i 

and  characteristic  function 

fx(?)  =  (l-i2T)*K.  (36) 

Since  random  variable  x  is  obviously  nonnegative  by  (35),  we  can  choose  shift 
b=0;  i.e.  y=x.  A  plot  of  the  cumulative  and  exceedance  distribution  functions 
of  random  variable  y  obtained  from  characteristic  function  (36)  with  K=4  is 
given  in  figure  3  for  0  <  v  £  2»/a.  The  values  of  a  and  L  have  been  chosen 
such  that  aliasing  and  tail  error  are  insignificant. 

The  ordinate  scale  for  figure  3  is  a  logarithmic  one.  The  lower  right 
end  of  the  exceedance  distribution  function  curve  decreases  smoothly  to  the 
region  1E-U,  where  round-off  noise  is  encountered.  The  exceedance 
distribution  function  values  continue  to  decrease  with  v  until,  finally, 
negative  values  (due  to  round-off  noise)  are  generated.  For  negative 
probability  values,  the  logarithm  of  the  absolute  value  is  plotted,  but 
mirrored  below  the  IE-12  level.  These  values  have  no  physical  significance, 
of  course;  they  are  plotted  to  illustrate  the  level  of  accuracy  attainable  by 
this  procedure  with  appropriate  choices  of  a  and  L. 

For  this  example,  Nsl/a=2666,  while  M=256.  Thus  collapsing,  according  to 
(21)  or  (32),  by  over  a  factor  of  10  has  been  employed  and  a  small  size  FFT 
has  been  utilized.  Nevertheless  the  error  realized  for  the  cumulative  and 
exceedance  distribution  functions  is  in  the  IE-12  range,  the  limit  of  accuracy 
of  the  Hewlett  Packard  9845B  Desk  Calculator  used  here.  Finer  spacing  in  the 
distribution  outputs  is  achievable  by  merely  increasing  M. 
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2.  Gaussian 


The  characteristic  function  for  a  zero-mean  unit-variance  random  variable 
is 


fx(j)  =  exp(-J2/2),  (37) 

and  the  probability  density  function  and  cumulative  distribution  function  are 
(ref.  11,  eq.  10.5.3) 


Px(v)  =  (2*)-172  exp(-v2/2),  Px(v)  -f(v).  (38) 


For  b  =  5ir/2,  using  (4), 


Py(0)  =  Px(-b)  =  <§(-b)  .  2E-15.  (39) 

which  is  negligible,  as  desired. 

Plots  of  the  cumulative  and  exceedance  distribution  functions  for  random 
variable  y  are  given  in  figure  4  for  1=7,  a=.3.  The  logarithmic  ordinate 
gives  rise  to  the  characteristic  parabolic  shape  on  the  tails  of  the 
distributions.  Once  again,  the  probabilities  decrease  to  the  level  of  the 
round-off  noise  and  fluctuate  around  IE— 12  near  the  edges  of  the  fundamental 
aliased  interval  (0,2Wa).  The  fact  that  the  cumulative  distribution  function 
of  y  starts  in  the  round-off  noise  at  v=0  indicates  that  b=5ir/2  was  large 
enough  to  guarantee  y  >  0  with  probability  virtually  1,  Also  indicated  on  the 
figure  is  the  origin  for  random  variable  x.  We  have,  from  (4), 

Px(u)  »  Py(u+b);  (40) 

thus  for  example 

Prob(x  <  0)  «  Px(0)  =  Py(b)  =  .5.  (41) 


14 


f 

4.  Gaussian;  L»7,  A". 3,  b»2.5Tr,  M*258 
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In  figure  5,  the  only  change  is  to  decrease  limit  L  from  7  to  6.  The 
tail  error  mentioned  in  (28)  et  seq.  then  dominates  the  round-off  noise  and 
lias  a  sinusoidal  variation.  Aliasing  is  not  a  problem,  as  witnessed  by  the 
fact  that  the  cumulative  and  exceedance  distribution  functions  of  random 
variable  y  have  decayed  below  IE-12  well  before  the  edges  of  the  interval  are 
reached. 

When  limit  L  is  restored  to  7,  and  sampling  increment  a  is  increased  to 
.5,  aliasing  becomes  significant,  as  shown  in  figure  6.  The  exceedance 
distribution  function  has  not  yet  decayed  to  the  round-off  noise  level  at 
v=2ir/a,  and  the  cumulative  distribution  function  shows  a  large  negative 
probability  region  near  v=0.  Shift  b  has  been  maintained  at  the  value  5w/2, 
corresponding  to  (39). 

When  L  and  a  are  restored  to  their  values  7  and  .3  as  for  figure  4,  but  b 
is  decreased  to  5«/3,  the  probability  of  y  becoming  negative  is,  from  (4)  and 
(38).  J  (-5»/3)  =  .82E-7.  This  is  reflected  in  the  cumulative  distribution 
function  for  y  in  figure  7  at  v*=G,  where  the  probability  value  is  well  above 
the  round-off  noise  level.  Also,  the  exceedance  distribution  function 
develops  significantly  negative  values  near  v  >»  2*/a. 

Accurate  evaluation  of  the  cumulative  and  exceedance  distribution 
functions  can  only  be  achieved  when  L,  A,  and  b  are  properly  chosen.  Probably 
the  optimum  combination  for  the  Gaussian  variate  is  displayed  in  figure  8, 
wt.ere  a  has  been  increased  to  .4,  the  distributions  are  centered  on  the 
fundamental  aliased  interval  (0,  2s/a)  by  choice  of  b,  and  l  is  taken  at  7  to 
avoid  tail  error. 

3.  Smirnov 

The  limiting  characteristic  function  of  a  measure  of  goodness  of  fit 
based  on  the  sample  distribution  function  was  derived  by  Smirnov  and  is  given 
by  (ref.  12,  eq.  30.104) 

V?)  -  (jr^iy)1/2  where  s  *  for  Ji °-  («) 
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An  expansion  about  J=0  yields 

fx($)  3  l+i  i.e.,  ux  =  1/6,  a3  =  1/45.  (43) 

And  since  the  goodness  of  fit  is  always  positive,  random  variable  x  is 
positive  and  we  can  choose 


b=Q . 


(44) 


Since 


sin((l+i)/?)~  \j  exp(VF(l-i))  as  £♦+*,  (45) 

it  follows  that 

fx(?)~  23/4f1/4  exp(-iyf+  U^VF-f))  as?-*-*.  (46) 

The  phase  of  this  term  rotates  according  to  V?1/ 2;  if  we  were  to  choose  b^O, 
fyllj)  would  rotate  faster  than  fx(?)  (linear  with  J  rather  than^f1).  This 
could  necessitate  a  faster  sampling  rate,  which  is  undesirable. 

The  cumulative  and  exceedance  distribution  functions  are  plotted  in 
figure  9.  L  and  a  have  been  chosen  so  as  to  avoid  tail  error  and  aliasing. 

The  exceedance  distribution  function  is  seen  to  decay  exponentially  until  it 
reaches  approximately  2E-11;  the  bump  in  the  curve  at  this  point  is  a 
manifestation  of  the  limited  accuracy  of  the  trigonometric  functions  built 
into  the  calculator  employed.  Larger  values  of  v  lead  to  round-off  noise 
around  the  IE— 12  level. 

A  comparison  of  results  for  this  characteristic  function,  with  Requicha's 
method  described  in  (29)  et  seq.,  is  given  in  figure  10  for  FFT  size  H-1Q24. 
The  plot  labeled  with  N*l=512  is  precisely  Requicha's  method.  Aliasing  is 
known  to  be  insignificant  for  4=1,  as  seen  by  reference  to  figure  9  and 
observing  that  extrapolation  of  the  straight  line  section  of  the  exceedance 
distribution  function  would  result  in  probability  values  near  IE— 13  at 
v=2 »/a.  The  dashed  portion  of  the  N=l=512  curve  in  figure  10  in  fact 


21 


TR  No.  7023 


corresponds  to  negative  probability  estimates;  these  grossly  inaccurate 
results  are  due  to  an  inadequate  value  of  limit  L,  leading  to  large  tail  error. 

When  N  is  simply  increased  to  1023,  the  middle  curve  in  figure  10  results 
from  Requicha's  method.  Again,  negative  estimates  are  indicated  by  the  dashed 
portion  of  this  curve,  although  two  orders  of  magnitude  smaller  than  above. 

The  reasons  for  these  errors  have  been  delineated  in  (29)  et  seq. 

The  bottom-most  curve  in  figure  10  (solid  curve)  is  that  obtained  by  the 
method  proposed  in  this  report  for  L  =  1023.  Exceedance  distribution  function 
estimates  in  the  IE— 10  range  are  obtained,  but  the  error  returns  to  the  IE— 8 
range  at  v=2ir/a.  No  negative  probability  values  occur.  Also,  by  simply 
increasing  limit  L,  while  keeping  FFT  size  M  fixed,  the  error  can  be  reduced 
significantly  further,  as  already  witnessed  by  figure  9. 

4.  Noncentral  Chi-Square 

Here  the  random  variable  x  is  given  by 

x  -  Jl  (gk  +  dk)2  (47) 

k=l 


where  [d^  are  constants,  and  {g^]  are  independent  Gaussian  random 

variables  with  zero-mean  and  unit  variance.  The  characteristic  function  of  x 

is 


fx(5)  -  (l-i2*rK/2  exp^f)  ,  (48) 

where  deflect  i-*n  d  is  defined  according  to 

d2  =  ^  dk’  (49) 

k-1 

We  actually  consider  a  more  general  characteristic  function  than  (48),  namely 
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fx(?)  =  (l-i2J)"v  exp^!^  =  exp(plr  '  vjtn(l-i2jr)),  (50) 

where  v  is  an  arbitrary  positive  real  constant.  Suppose  that  we  use  the 
principal  value  logarithm  for  J?n(z),  where  the  branch  cut  lies  along  the 
negative  real  axis  of  the  complex  z  plane  (ref.  13,  sect.  4.1.1).  Then  since 
the  argument  of  the  logarithm  in  (50)  never  crosses  the  branch  cut,  form  (50) 
gives  the  correct  characteristic  function  values  automatically  for  all  real  9, 
and  any  v. 

The  probability  density  function  and  exceedance  distribution  function 
corresponding  to  (50)  are  (ref.  14,  6.631  4) 

px(v)  1  Iv-l<d1'7)  for  v>0’ 

4*40 

1  -  Vv)  -  Idt  t  exp f-  1  Iv_t(dt)  *  Qv(d,V\T)  for  v  >  0.  (51) 

Since  the  probability  density  function  in  (51)  is  never  negative  (ref.  13, 
sect.  9.6.1),  (50)  is  a  legal  characteristic  function.  Also  because  random 
variable  x  is  always  positive  according  to  (51),  we  choose  shift  b=0.  Plots 
of  the  exceedance  distribution  function,  as  determined  from  characteristic 
function  (50)  are  displayed  for  various  values  of  d  in  figure  11.  The  values 
of  L  were  chosen  for  each  d  value  so  as  to  control  the  tail  error  below  the 
IE— 10  level  plotted.  Direct  calculation  of  the  exceedance  distribution 
function  directly  from  (51)  would  be  a  formidable  task  for  arbitrary  v  values. 

5.  Product  of  Correlated  Gaussian  Variates 

Let 

X  a  st  (52) 

where  s  and  t  are  zero-mean  unit-variance  Gaussian  random  variables  with 
correlation  coefficient  p.  The  joint  probability  density  function  of  s  and  t 
is 
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Pct(u,v)  =  (2nYl-p  J  exp  - 


2  2 

u+v  -2puv 
2(1-pZ) 


The  characteristic  function  of  x  is  then 

fx(£)  =  exp(ifst)'  =  jjdu  dv  exp(ifuv)  pst(u,v)  = 

-  l-i2p$+(l-p2)f  2j  U2  =  [l-1(l+p )f]~1/2  [l+i(l-p)*]"1/2,  (54) 

via  repeated  use  of  ref.  14,  eq.  3.323  2.  The  corresponding  probability 
density  function  of  x  is 

p*(v)  "  '^ex{€?)  K°(&)  for  a"  v'  <55) 
via  ref.  14,  eq.  3.478  4. 

(If  we  transform  this  probability  density  function  according  to  (1)  and 
use  ref.  14,  eq.  6.611  9  and  ref.  13,  eq.  4.4.15,  we  get  precisely  (54). 
Alternatively,  if  we  transform  (54)  and  modify  the  contour  to  wrap  around  the 
branch  line  along  the  imaginary  axis  and  then  use  ref.  14,  eq.  3.388  2,  we  get 
(55).  Or  we  can  use  ref.  14,  eq.  3.754  2.) 

We  actually  consider  a  more  general  characteristic  function  than  (54), 
namely 


{(f)  «  [l-i2prU-o2)52']  V  -  exp^n[l-i2p 


The  mean  of  this  random  variable  x  is  given  by 


J  +  (l-0C) 


a  2vp. 
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The  probability  density  function  corresponding  to  (56)  is 


px(v)  =  ^dj  exp(-i?v)  fx(?)  = 


We  can  move  the  contour  in  (58)  to  the  real  y-axis,  because  the  branch  points 
of  the  integrand  are  at  y  =>  ±  which  are  outside  the  path  of 

integration,  since  ]o \  <  1.  Then  using  ref.  14,  eq.  3.771  2  and  ref.  13, 
eq.  6.1.17,  we  obtain 


Since  this  probability  density  function  is  never  negative  (ref.  13,  sect. 
9.6.1),  (56)  is  a  legal  characteristic  function.  If  we  Fourier  transform  (60) 
via  ref.  14,  6.699  12,  we  get  (56)  directly. 


There  is  no  simple  relation  for  the  cumulative  distribution  function  of 
this  random  variable.  Nevertheless,  it  is  a  simple  matter  to  evaluate 
directly  from  characteristic  function  (56).  The  Jtn  in  (56)  causes  no  problems 
since  its  argument  never  crosses  the  branch  cut.  A  plot  for  \>«7.7  and  o»-.3 
is  displayed  in  figure  12.  The  rate  of  decay  of  the  distribution  is  different 
for  each  tail.  The  round-off  noise  is  clearly  visible  at  both  ends  of  the 
range  of  v  values. 
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We  now  have  the  capability  to  handle  the  following  type  of  statistical 
problem  in  a  fairly  easy  fashion.  Consider  random  variable 


x 


(61) 


where  £r^  are  arbitrary  random  variables,  statistically  independent  of  each 
other,  and  with  different  distributions.  Power  \>k  is  arbitrary  (except  that 
must  be  a  positive  integer  for  those  rk  that  can  become  negative).  Let 
the  probability  density  function  of  random  variable  rk  be 

characteristic  function  of  r^  is 


pk(v).  Then  the 


gk(?)  «  exP^rkk)  3  Jdv  exp(i?v  k)  P|JV)  * 
-  exp(,r>  >U'’k  Pi<(t1,''k)  • 


(62) 


If  (62)  is  rot  integrable  in  closed  form,  it  can  be  evaluated  by  means  of  an 
FFT  (one  for  each  k  if  the  probability  density  functions  or  vk  are  all 
different).  Then  toe  characteristic  function  of  random  variable  x  in  (61)  is 
given  by 


K 


f*if>  -tm  • 

k*l 


(63) 


Now  the  techniques  of  this  report  are  directly  applicable  to  (63). 
An  additional  example  is  afforded  by 
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ak  vk 


*  2h\  '2* 


where  {a^,  {0^,  and  (yk^  are  constants,  and  [vk^  are  independent 
random  variables  with  arbitrary  probability  density  functions.  The 


characteristic  function  of  x  s 


fx(f)  =  exp(iTx)  =  JdV  pv(V)  exp(if£akvk  +  £skvk)2  +^ykvk]),  (65) 


where  V  =  (v^  vK).  Now  sir 


(~j  Jdy  exp(-iay2  ♦  iby)  =.  exP^fy  for  a  *  0,  (66) 

we  identify  a  =  fl 4,  b  Bkvk,  eliminate  the  square  in  the  exponent, 
and  express  (65)  as 

%(?)  «  J*dV  Pv(V)  exp(ij^akv2  ♦  if^vO* 

*(t£)  {  dy  exp^  ♦  ^v^Vk)  - 

5  K 

"  ($)  Jdy  exp(~  ^~)  [|dvk  pk(V  exp(i*(Vk  *  Vk  *  (67) 


where 


pv(,)  ■Tki'li)  •  <«: 

k-1 

The  inner  integrals  in  (67)  can  either  be  done  analytically  or  numerically. 
Then  the  remaining  single  integral  on  y  must  be  numerically  evaluated  to  find 
characteristic  function  fx(f).  As  an  example,  if  vk  is  exponentially 
distributed 
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Pk(v)  =  ak  exp(-akv)  for  v  >  0, 

then  the  inner  integrals  in  (67)  are  w-functions;  see  ref.  13,  ch.  7.  A 
simpler  method  of  handling  general  quadratic  expressions  like  (64)  with 
Gaussian  V  is  presented  in  ref.  15. 
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SUMMARY 

An  accurate  method  for  efficient  evaluation  of  the  cumulative  and 
exceedance  distribution  functions  has  been  derived  and  applied  to  several 
examples  to  illustrate  its  utility.  Choice  of  the  sampling  increment  a 
applied  to  the  characteristic  function  controls  the  aliasing  problem,  and 
selection  of  the  limit  L  minimizes  the  tail  error;  the  effects  of  both  of 
these  parameters  can  be  observed  from  sample  plots  of  the  distributions  and 
can  be  modified  if  needed.  Additionally,  shift  b  must  be  chosen  so  as  to 
yield  a  positive  random  variable  with  probability  virtually  1.  The  number  of 
distribution  values  yielded  depends  on  the  size  of  the  FFT  employed  and  can  be 
independently  selected  to  yield  the  desired  spacing  in  distribution  values. 
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APPENDIX  A.  SAMPLING  FOR  A  FOURIER  TRANSFORM 
Suppose  we  are  interested  in  evaluating  Fourier  transform 
G(f)  =  jdt  exp(-i2irft)  g(t). 


(A-l) 


If  we  sample  at  interval  a  in  t  in  (A-l),  and  use  integration  weighting  w(t), 
we  have  the  approximation  to  G(f), 


0(f)  h 


dt  exp(-i2*ft)  g(t)  S  (t)  w(t) 


=  G(f}<&  K(f) 


=  l2G(f~ 


(A-2) 


where  infinite  impulse  train  (sampling  function) 


SA(t)  = 


(A-3) 


and  ® denotes  convolution. 


The  term 


HG(f-  <A-4) 

n 

in  (A-2)  is  an  infinitely  aliased  version  of  desired  function  G(f);  this 
aliasing  is  an  unavoidable  effect  due  to  sampling  at  increment  a.  However,  to 
minimize  any  further  aliasing  in  (A-2),  we  would  like  W(f)  =  £(f),  which 
requires  w(t)  =  3.  for  all  t;  strictly,  all  we  need  is 

w(na)  3  1  for  all  n.  (A-5) 


A-l 


...  .  ••y.v-A 


•  -v  ri  *  : 
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That  is,  the  best  weighting  in  (A-2)  is  uniform. 


As  an  example,  for  Simpson's  rule,  we  have  weighting 

w(na)  =  ...,  -j,  y,  j,  -j,  ...  =  1  +  -j(-l)n  or  1  -  ■j(-l)n, 

which  can  be  represented  as  samples  of  time  function 

w(t)  »  1  +  exp ( i -rrt/ a )  or  1  -  exp(iirt/a). 


The  corresponding  transform  is 


w(f|  ■  I 


dt  exp( -i 2-rrf t)  w(t)  = 


=  $(f)  +  -J  "  7^)  or  a(f)  -  -J  S(f  ~  ^-)  • 


(A-6) 


(A- 7 ) 


( A— 8 ) 


But  this  window  function  substituted  in  (A-2)  results  in  an  extra  aliasing 
lobe  in  G(f),  halfway  between  the  unavoidable  major  lobes  of  (A-4)  at 
multiples  of  i/a,  of  magnitude  1/3  as  large.  This  effect  very  adversely 
affects  the  quality  of  6(f)  insofar  as  its  approximation  to  the  desired  G(f) 
is  concerned.  Thus  the  best  sampling  plan  in  (A-2)  is  the  equal  weight 
structure  of  (A— 5 )  when  one  wants  to  approximate  the  Fourier  transform  of 
(A-l).  For  a  bounded  region,  this  is  modified  to  the  Trapezoidal  rule,  i.e., 
half-size  weights  at  the  boundaries. 
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APPENOIX  B.  LISTINGS  OF  PROGRAMS  FOR  FIVE  EXAMPLES 

The  following  listings  are  programs  in  BASIC  for  the  Hewlett  Packard 
9845B  Desktop  Calculator.  The  FFT  utilized  is  one  with  the  capability  of  a 
zero  subscript  and  is  listed  at  the  end  of  the  appendix.  Mathematically,  the 
FFT  programmed  is 


M-l 

Zffl  »  ^  exp(-i2irmk/M)  z^  for  0  £  m  <  M-l, 
k=0 

r  nM-1  c  iM-1 

where  the  arrays  and  are  handled  directly,  including  the  zero- 

subscript  terms  zQ  and  ZQ. 

A  detailed  explanation  of  the  first  program  below  for  Chi-Squared  random 
variables  is  as  follows:  line  20  specifies  the  parameter  K,  where  2K  is  the 
number  of  squared-Gauss i an  random  variables  summed  to  yield  random  variable 
x.  Lines  30-60  require  inputs  L,  a,  b,  M  respectively,  on  the  part  of  the 

user.  Line  110  is  the  input  of  mean  ux  of  random  variable  x,  as  evaluated 

analytically  from  characteristic  function  fx($).  Lines  180-210  specifically 
evaluate  the  characteristic  function  fy(£)  at  general  point?.  All  of  these 
lines  mentioned  thus  far  require  inputs  on  the  part  of  the  user  and  are  so 

noted  in  the  listing  by  the  presence  of  a  single  I  on  each  line;  the  comments 

after  a  douDle  I!  are  for  information  purposes  only  and  need  not  be  modified. 
This  convention  is  also  adopted  in  the  remaining  listings. 

Lines  220-240  accomplish  the  collapsing  operation  of  ( 32 ) — ( 33) .  The 
cumulative  and  exceedance  distribution  functions  are  finally  evaluated  and 
stored  in  arrays  X(*)  and  Y(*)  in  lines  400-410. 

Some  further  elaboration  is  necessary  for  the  listing  of  the  Smirnov 
characteristic  function  as  given  by  (42).  Since  a  characteristic  function  is 
a  continuous  function  of  real  f ,  the  square  root  in  (42)  is  not  a  principal 
value  square  root,  but  in  fact  must  yield  a  continuous  function  in  f.  In 
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order  to  achive  this,  the  argument  of  the  square  root  is  traced  continuously 
from  $ =0  (line  110).  If  an  abrupt  change  in  phase  is  detected,  a  polarity 
indicator  takes  note  of  this  fact  (line  250)  and  corrects  the  final  values  of 
characteristic  function  fy(f)  (lines  260-270).  No  problems  are  encountered 
with  complex  sin(z)  since  it  is  analytic  for  all  z. 


18  ! 
20 
30 

40 

S0 

60 

70 

80 

90 

100 

110 

123 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 


CHI -SQUARE 
K=4 
L=*200 

Del  t.  i=,  075 

Bs=fi 

M=2-'-8 


CHARACTERISTIC  FUNCTION  l/tl-i  2  xi )M 
!  2K=8  degrees  of  freedom 

!  Limit  on  integral  of  char. 

!  Sampling  increment  on  char. 
!  Shift  b 
!  Size  of  FFT 


f unc t i on 
funct i on 


PRINTER  IS  0 

PRINT  "L  = ” ; L , "Delta  ; Del ta, “b 
REDIM  Xf0:M-l),Yt0:M-l) 

DIM  X<0: 1023) ,  Yt0: 1023) 

Mux=2*K 
Muy=Mux+Bs 
X  <  8  )  =>0 

Y<0>». 5*De 1 t a*Muy 
N= I  NT  <  L/De 1 1  a) 

FOR  Ns=>  1  TO  N 
Xi =Del t  a*Ns 
C=X i +X  i 

CALL  Multi, -C, 1 , -C, A, B) 

CALL  Mul <A,B,A,B,C,D) 

CALL  Di  v( 1 , 0, C , D, Fyr , Fy 1  ) 

Ms=Ns  MOD  M 
X  <  Ms )»X  <  Ms >+Fyr> Ns 
Y(Ms)=YtMs)+Fyi  •'Ns 
NEXT  Ns 

CALL  Fftl0z<M.Xt*>,Y<*)> 


=  “ J  Bs, “M  ="JM 


!  Mean  of  random  variable  x 


Argument  xi  of  char.  fn. 
Calculation  of 
charac  t er i st i c 
f unct ion  fytxi  > 
for  K=4 
Co  1 1 apsi ng 


1  !  0  subscr i pt  FFT 
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270  PLOTTER  IS  "GRAPHICS" 

280  GRAPHICS 

290  SCALE  0,11,-14,0 

300  LIME  TYPE  3 

310  GRID  M/8,  1 

320  PENUP 

330  LINE  TYPE  1 

340  B=Bs*M*Del ta/<2*PI > 

350  MOVE  B , 0 

360  DRAW  B, -14 

370  PENUP 

380  FOR  Ks=0  TO  M-l 

390  T=Y<Ks)/Pl-Ks/M 

400  X<Ks>*. 5-T 

410  Y<Ks>=Pr=. 5+T 

420  IF  Pr>=lE-12  THEN  Y=LGTCPr> 

430  IF  PrO-lE~12  THEN  Y=-24-LGT<-Pr > 

443  IF  ABS<Pr)< IE-12  THEN  Y  =  -12 

450  PLOT  Ks,Y 

460  NEXT  Ks 

470  PENUP 

480  PRINT  Y<0>5Y<1>;Y<M-2>;Y<M-1> 

490  FOR  Ks*0  TO  M-l 

503  Pr*X<Ks ) 

510  IF  Pr>=lE-12  THEN  Y=LGT ( Pr > 

520  IF  Pr<*-1E~12  THEN  Y=-24-LGT < -Pr ) 

530  IF  ABSCPrX  IE-12  THEN  Y«-12 

540  PLOT  Ks , Y 

550  NEXT  Ks 

563  PENUP 

570  PAUSE 

580  DUMP  GRAPHICS 

590  PRINT  LINC5) 

600  PRINTER  IS  16 

610  END 

620  ! 

630  SUB  MuHXl.Yl,X2,Y2,A,B> 

640  A*X 1 *X2-Y 1 *Y2 

650  B«X1*Y2+X2*Y1 

660  SUBEND 

670  ! 

680  SUB  D  i  v  ( X 1 ,  Y 1 ,  X2 ,  Y2 ,  A ,  B  > 

690  T  ®X2*X2+Y2*Y2 

700  A«<.X1*X2+Y1*Y2->'T 

710  B*<Y1*X2-X1*Y2>T 

720  SUBEND 

730  ! 

740  SUB  Ff t lOz < N, X ; * > , Y< * ) »  <  N  <= 


Origin  for  random  variable  x 


Cumulative  probability  in  X<*> 
Exceedance  probability  in  Y(*> 


Z  1*22 


21  ■'22 


-'•13  a  1024,  N=2A I NTEGER  0  subscript 
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10  !  GAUSS  I AN  CHARACTERISTIC  FUNCTION  exp<-.5  xiA2) 

20  L  =  7  !  Limit,  on  integral  of  char,  function 

30  Delta®. 3  !  Sampling  increment  on  char,  function 

40  Bs=. 375*<2*PI/Del ta)  !  Shift  b,  as  fraction  of  alias  interval 

50  11=2^3  !  S  i  ze  of  FFT 


60  PRINTER  IS  0 

70  PRINT  "L  =";L, "Delta  =";Delta,"b 
80  REDIM  X<0:M-1), Y<0:M-1) 

90  DIM  X<0: 1023), Y<0: 1023) 

100  Mux=0 

110  Muy=Mux+Bs 

120  X<0)=0 

130  Y<0)=. 5*Del ta*Muy 

140  N=INT<L/Delta) 

150  FOR  Ns=l  TO  N 

160  X i =De 1 1  a*Ns 

170  A=EXP  < - . 5*X i  *X i  ) 

180  B=Bs*Xi 

190  Fyr=A*C0S<B) 

200  Fy  i ®A*SIN<B) 

210  Ms=Ns  MOD  M 

220  X<Ms)=X<Ms)+Fyr/Ns 

230  Y<Ms)=Y<Ms)+Fyi/'Ns 

240  NEXT  Ns 

250  CALL  Fftl0Z<M,X<*), Y<*>) 

260  PLOTTER  IS  "GRAPHICS" 

270  GRAPHICS 

280  SCALE  0, M, - 1 4, 0 

290  LINE  TYPE  3 

300  GRID  M/8,  1 

310  PENUP 

320  LINE  TYPE  1 

330  B=Bs*M*Del ta/(2*PI > 

340  MOVE  B ,  0 

350  DRAW  B , - 1 4 

360  PENUP 

370  FOR  Ks=0  TO  M-l 

380  T=Y<Ks)/PI-Ks/M 

390  X<Ks ) * . 5-T 

400  Y<  Ks ) ®Pr® . 5+T 

410  IF  Pr  >* 1 E- 1 2  THEN  Y=LGT <Pr> 

420  IF  Pr <  =- 1 E- 1 2  THEN  Y=-24-LGT < -Pr ) 

430  IF  ABS < Pr )< 1 E- 1 2  THEN  Y=-12 

440  PLOT  Ks ,  Y 

450  NEXT  Ks 

460  PENUP 

470  print  y<0);yu>;y<m-2);y<m-i) 

480  FOR  Ksa0  TO  M-l 

490  Pr®X(Ks) 

500  IF  Pr >*1E-12  THEN  Y»LGT<Pr) 

510  IF  Pr  <  «*-  1  E-  t  2  THEN  Y*-24-LGT v -Pr > 

520  IF  ABS< Pr ) < 1 E  —  1 2  THEN  Y®-12 

530  PLOT  Ks,  Y 

540  NEXT  Ks 

550  PENUP 

560  PAUSE 

570  DUMP  GRAPHICS 

580  PRINT  L I N < 5 ) 

590  PRINTER  IS  16 

600  END 

610  ! 

620  SUB  F  f  1 1  0  z  C  N ,  X  •.  *  > ,  Y  >  *  >  *  !  N  -  = 

B-4 


=  " ;  Bs, 11  m  =";m 


!  Mean  of  random  variable  x 


Argument  xi  of  char.  fn. 
Calculation  of 
character  i  st  i  c 
f unc  t i on 
fy(xi ) 

Col  1  aps i ng 


1  !  0  subsc  r i pt  FFT 


M  Origin  for  random  variable  x 


1!  Cumulative  probability  in  X < * ) 
!!  Exceedance  probability  in  Y<*) 


2' 1 0  *  1024,  N»2MNTEGER 


0  subscript 
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10  !  SMIRNOV  CHARACTERISTIC  FUNCTION  C s/s i n < % ) ] A 1 /2  where  s=(  1  +  i )sqr <xi ) 


20 

L=3@00  !  Limit  or 

integral  of  char,  function 

30 

D  e  1 1  a= 1  !  Sampl i ng 

incr-.ment  on  char,  function 

40 

Bs=0  !  Shift  b 

50 

M=2A8  !  Size  of  FFT 

€0 

PRINTER  IS  0 

70 

PRINT  "L  =";L,  “Delta  ='* ;  Del  ta,  "b 

_  u 
—  J 

as, "M  =";m 

80 

REDIM  X<0:M-1),Y(0:M-1) 

90 

DIM  X<0! 1023), Y<0! 1023) 

100 

Mux=l/6 

| 

Mean  of  random  variable 

X 

110 

R=0 

1 

Argument  of  square  root 

120 

P=1 

1 

Polarity  indicator 

130 

Muy=Mux+Bs 

140 

X  <  0  )  3  0 

150 

Y(0)=.5*Delta*Muy 

160 

N=INT<L/De1ta) 

170 

FOR  Ns=l  TO  N 

180 

Xi =Del ta*Ns 

j  i 

Argument  xi  of  char,  fn 

. 

190 

A=SQR  <  X i ) 

j 

Cal cul at i on 

200 

CALL  Sin<A,A,B,C) 

j 

of 

210 

CALL  Di v(A, A, B, C, D, E) 

character i st i c 

220 

CALL  Sqr<D, E, A, B) 

j 

f  unc  t i on 

230 

Ro=R 

! 

fyCxi  ) 

240 

R=ATN(B/A) 

! 

250 

IF  ABS < R-Ro ) > 1 . 6  THEN  P=-P 

! 

260 

Fyr=A*P 

i 

270 

Fy i =  B*P 

i 

280 

Ms=Ns  MOD  M 

t  t 

Col  1 apsi ng 

290 

X<Ms)=X<Ms)+Fyr/Ns 

300 

Y<Ms)=Y<Ms)+Fyi/Ns 

310 

NEXT  Ns 

320 

CALL  Fftl0z<M,X<*),  ■{<*>  ) 

i  | 

0  subscript  FFT 

330 

PLOTTER  IS  "GRAPHICS" 

340 

GRAPHICS 

350 

SCALE  0 , H , - 1 4 , 0 

360 

LINE  TYPE  3 

370 

GRID  M/8, 1 

380 

PENUP 

390 

LINE  TYPE  1 

400 

B=Bs*M*Del t  a-  < 2*P I > 

i  i 

Origin  for  random  variable 

X 

410 

MOVE  B , 0 

420 

DRAW  B,-U 

430 

PENUP 

440 

FOR  ks=»0  TO  M- 1 

450 

T«Y<Ks)/PI-Ks/M 

460 

;:<Ks)«,  5-T 

t  t 

Cumulative  probability 

\  n 

X<*) 

470 

Y<Ks)=»Pr  =  .  5+T 

i  i 

Exceedance  probability 

i  n 

Y<  *  > 

480 

IF  Pr  >■ 1 E- 1 2  THEN  Y  =  LGT < Pr > 

490 

IF  Pr<  —  IE-12  THEN  Y=«-24-LGT< -Pr  > 

500 

IF  ABSCPrls IE-12  THEN  Y»-12 

510 

PLOT  Ks ,  Y 

520 

NEXT  Ks 
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530  PENUP 

540  PRINT  Y<0> ; YC 1  > ; Y<M-2>  J  YCM-l > 

550  FOR  Ks=0  TO  M-l 

560  Pr=X<Ks> 

570  IF  Pr>=  IE- 12  THEN  Y  =  LGT <Pr > 

530  IF  Pr<=-1E-12  THEN  Y=-24-LGT(-Pr > 

590  IF  RBS<PrXlE-12  THEN  Y=-12 

600  PLOT  Ks , Y 

610  NEXT  Ks 

620  PENUP 

630  PAUSE 

640  DUMP  GRAPHICS 

650  PRINT  LIN<5) 

660  PRINTER  IS  16 

670  END 

680  ! 

690  SUB  Dw(Xl,Yl,X2,Y2,A,B)  !  Z1/Z2 

700  T»X2*X2+Y2*Y2 

710  A=<X1*X2+Y1*Y2XT 

720  B=<Y1*X2-X1*Y2)'/T 

730  SUBEND 

740  ! 

750  SUB  Sqr (X, Y , A , B)  !  PRINCIPAL  SQR(Z) 

760  IF  XO0  THEN  800 

770  A=B=SQRt.5*RBS<Y>) 

780  IF  Y< 0  THEN  B=-B 

790  GOTO  910 

800  F=$QR<$QR<X*X+Y*Y)> 

810  T=,  5*ATN<  Y-'X) 

820  A=F*COS<T) 

830  B=>F*SIN<T> 

840  IF  X>0  THEN  910 

850  T=R 

860  A*-B 

870  B=*T 

880  IF  Y>«0  THEN  910 

890  A=-A 

900  B  — B 

910  SUBEND 

920  ! 

930  SUB  StruX,  Y,  A,  B>  •  SIN(Z) 

940  E*EXP<Y) 

950  A-.5*SIN<X)*(E+XE) 

960  IF  ABS< Y)< . 1  THEN  990 

970  S».5*<E-1/E) 

980  GOTO  1010 

990  S*Y*Y 

1000  S-Y*U20*S*^20-*-S‘>  120 

1010  B«COS<X  >*$ 

1020  SUBEND 

1030  i 

1O40  SUB  Fft  10:<N,  ;<(»),  Vv«'  )  i  N  2X0  *  1024,  N-2MHTEGER  0  subscript 


B-6 


TR  No.  7023 


NON-CENTRAL  CH I -SQUARE  CHARACTERISTIC  FUNCTION 
exp<i  d'-2  xi  x  s>  x  s*nu  where  s  =  1  - i  2  xi 


Nu*2.7  !  Power  law  nu 

Ds=3  !  Deflection  d 

L=500  !  Limit  on  integral  of  char,  function 

Delta=.05  !  Sampling  increment  on  char,  function 

Bs=0  !  Shift  b 

M=2/'8  !  Size  of  FFT 

PRINTER  IS  0 

PRINT  "L  =";L,  "Delta  =*" ;  Del  ta,  "b  =“;Bs,"M  =  ";M 
REDIM  X<0:M-1),Y<0:M-1> 

DIM  X<0: 1023), Y(0: 1023) 

D2=Ds*Ds  !  Calculate 

Mux*2*Nu+D2  !  Mean  of  r 

Muy=Mux+Bs 
X<0)=0 

Y<0)=.5*Delta*Muy 
N*INT(L/Del ta) 

FOR  Ns»l  TO  N 

Xi=Delta*Ns  !!  Argument 

T=Xi +Xi  !  Cal cu 1  at i 

CALL  D i u<0, D2*Xi , 1 , -T, A, B)  !  character 

CALL  Log< 1 , -T, C, D)  !  function 

CALL  Exp<A-Nu*C, B-Nu*D+Bs*Xi , Fyr , Fyi )  !  fy(x 


Calculate  parameter 
Mean  of  random  variable  x 


Argument  xi  of  char.  fn. 
Cal cu 1  at i on  of 
character i st ic 
funct i on 

!  fy(xi) 


Ms=Ns  MOD  M 
X(Ms)=X<Ms)+Fyr/Ns 
Y(Ms)=Y(Ms)+FyixNs 
NEXT  Ns 

CALL  Fftl0z<M,X<*),Y<*)> 

PLOTTER  IS  “GRAPHICS" 

GRAPHICS 
SCALE  0 , M, - 1 4 , 0 
LINE  TYPE  3 
GRID  Mx8, 1 
PENUP 

LINE  TYPE  1 
B=Bs*M*D« 1 ta/ v  2*P  I > 

MOVE  B, 0 
DRAW  B, -14 
PENUP 

FOR  Ks*0  TO  M-l 
T*Y<Ks>/PI-Ks/M 
X(K*)-. 5-T 
Y’vKs  >«Pr«.  3+T 

IF  Pr  >* 1 E- 1 2  THEN  Y3LGT <  Pr  ' 

IF  Pr<  *- IE- 12  THEN  Y--24-LGT <. -Pr > 
IF  ABS <Pr)< IE-12  THEN  Y--12 
PLOT  Ks.Y 
NEXT  Ks 
PENUP 


Co  1 1 apsi ng 


8  subscript  FFT 


Origin  For  random  variable  x 


11  Cumulative  probability  in  X(*> 
11  Exceedance  probability  in  Yte> 
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319  PRINT  Y<0>;Y<1);Y<M-2>;Y<M-1> 

520  FOR  Ks*0  TO  M-l 

530  Pr*X<Ks> 

540  IF  Pr >*1E-12  THEN  Y-LGT<Pr> 

550  IF  Pr<  —  IE-12  THEN  Y--24-LGT < -Pr  ) 

560  IF  ABS<Pr)< IE-12  THEN  Y--12 

570  PLOT  Ks.Y 

580  NEXT  Ks 

590  PENUP 

600  PAUSE 

610  DUMP  GRAPHICS 

620  PRINT  LIN<5> 

630  PRINTER  IS  16 

640  END 

650  ! 

660  SUB  Diu<Xl,Yl,X2,Y2,A,B>  !  21/22 

670  T*X2*X2+Y2*Y2 

680  A»<X1*X2+Y1*Y2)/T 

690  B»<Y1*X2-X1*Y2>/T 

700  SUBEND 

710  ! 

720  SUB  Exp<X,Y,A,B>  !  EXP<2> 

730  T*EXP<X) 

740  A-T*C0S<Y) 

750  B*T*SIN< Y) 

760  SUBEND 

770  ! 

730  SUB  Log<  X, Y, A, B )  !  PRINCIPAL  LOG<Z> 

790  A».5*L0G<X*X+Y*Y> 

800  IF  X<  >0  THEN  830 

810  B".5*PI*SGN<Y) 

820  GOTO  830 

830  B«ATN<Y/X> 

840  IF  X<  0  THEN  B-B+P I *<t -2* < Y< 0 > > 

850  SUBEND 

860  ! 

870  SUB  Fftl0z<N,X<*>,Y<*>>  !  N  <*  2A10  ■  1024,  N*2A INTEGER 


O  subscript 
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10  • 

20  Nu*7. 

30  Rho*-.3 

40  L*5 

50  Delta*. 06 

60  Bs=.  5»<2*PI,'Del u) 

70  M*2^9 

80  PRINTER  IS  0 

90  PRINT  "L  *";L, "Delta 

100  REDIM  X<0:M-1),Y<0:M-1> 

110  DIM  X<0: 1023) , Y<0: 1023) 

120  T 1 = 1 -Rho*Rho 

130  T2*2*Rho 

140  Mux=2*Nu#Rho 

150  Muy*Mux+Bs 

160  X<0>*0 

170  Y<0)*.5*Delta*Muy 

130  N*INT  <LxDel ta) 

190  FOR  Ns* 1  TO  N 

200  Xi*De1ta*Ns 

210  CALL  LogU+Tt*Xi  *Xi  , -T2*Xi  ,  ft,  B) 

220  CALL  Exp<-Nu*A, Bs*Xi -Nu*B, Fyr , Fy i 

230  Ms*Ns  MOD  M 

240  X(Ms)*X<Ms)+FyrxNs 

250  Y<Ms)*Y<Ms)+Fyi/Ns 

260  NEXT  Ns 

270  CALL  Ff  t 10z<M,X<*),Y<*>) 

280  PLOTTER  IS  "GRAPHICS" 

290  GRAPHICS 

300  SCALE  0 , M , - 1 4 , 0 

310  LINE  TYPE  3 

320  GRID  Mx8, 1 

330  PENUP 

340  LINE  TYPE  1 

350  B*Bs*M*Delta-'<2*PI ) 

360  MOVE  B,0 

370  DRAW  B , - 1 4 

380  PENUP 


= " ; Bs , "M  *" ; M 


Cal cul ate 

parameters 

Mean  of  random  variable  x 


! !  Argument  xi  of  char.  fn. 

!  Calculation  of 

)!  character i st i c  function  fy(xi) 
!  !  Coll apsi ng 


!  !  0  subscript  FFT 


Origin  for  random  variable  x 


GAUSSIAN  PRODUCT  CHARACTERISTIC  FUNCTION  <56> 

7  !  Power  Nu 

Correlation  coefficient 
Limit  on  integral  of  char,  function 
Sampling  increment  on  char,  function 
Shift  b»  as  fraction  of  alias  interval 
Size  of  FFT 

"J Delta. “b 
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FOR  Ks«0  TO  M-l 
T-YCKsXPI-KsxM 
X<Ks>*. 5-T 
Y<k's>«Pr-.5+T 

IF  Pr>- IE- 12  THEM  >LGT<Pr) 

IF  Pr<»-1E-12  VHEn  Y«-24-LGT<-Pr> 
IF  ABS<Pr><  IE-12  THEN  Y— 12 
PLOT  K» , Y 
NEXT  K* 

PENUP 

print  y<0>;y<i>;y<m-2>;y<m-1) 

FOR  Ks*0  TO  M-l 
Pr=X<Ks> 

IF  Pr >- 1 E- 1 2  THEN  Y*LGT <Pr  > 

IF  Pr<--1E-12  THEN  Y--24-LGT< -Pr > 

IF  ABSCPrX  IE-12  THEN  Y—12 

PLOT  Ks,Y 

NEXT  Ks 

PENUP 

PAUSE 

DUMP  GRAPHICS 
PRINT  LIN<3> 

PRINTER  IS  16 
END 


H  Cumulative  probability  in  X<*> 
!!  Exceedance  probability  in  Y<*> 


SUB  Exp<X,Y,A,B> 
T*EXP<X) 

A*T*C0S< Y) 
B-T*SIN<Y) 

SUBEND 


!  EXP<Z> 


SUB  Log<X, Y, A, B> 
A-.3*L0G<X#X+Y*Y> 

IF  XO0  THEN  730 
B».3*PI*SCN<Y> 

GOTO  770 
B*ATN< Y^X) 

IF  X<0  THEN  B»B+PI*v 1-2*<Y<0>> 
SUBEND 


!  PRINCIPAL  L0G<Z > 


SUB  Fft lOz(N,X<#>,Y<*> 


•  N  <i 


1024,  N»2A INTEGER 


O  subscript 
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!  N  O  2 -s10  -  1024,  N«2"'  INTEGER 


0  subscript 


10  SUB  Fft  10z<N,X'*  V- 

20  DIM  C<0;236) 

30  INTEGER  II,  12,  13,  14,  15,  16, 17, 18, 19, 1 10, J,K 

40  DATA  1,  ,  99998 1 1 75283 , . 999924701339, . 999830581 796, . 9996933 1 8696, . 9995294 1 73 

01 , . 999322384588, . 999077727753, . 998793456205, . 998473580573, . 9931 13112900 

50  DATA  .997723066644, .997290456679, . 996320299291 ,. 996312612183, . 995767414468 

, .995184726672, . 994364570734, . 993906970002, . 993211949235, . 992479534599 

60  DATA  .991709733669, .990902633428, . 99005321 0262, . 9891 76509965, . 983237367731 

, .987301413153, .986308097245, .935277642389, . 98421 0092387,  .  983 1 0548743 1 

70  DATA  .9819638691 18, .980735280403, . 979569765685, . 97831 7370720, . 977028142658 

, .975702130039, . 974339382786, . 972939952206, . 971503890986, .  970031233195 

30  DATA  .963522094274, .966976471045, . 965394441698, . 963776065795, . 962121404269 

, .960430519416, .958703474896, .956940335732, . 955141 163306, . 933306040354 

90  DATA  .951435020969, . 949528 1 80593 ,. 94758559 1 01 3 , . 943607325381 ,. 943593458162 

, . 94 1544085 1 83, . 93°  "9223602, . 93733901 1913, . 935133509939, . 932992793835 

100  DATA  .9307669*1879, .923506080473, . 926210242138, .923379532511, .921514039342 

, .919113851690, .916679059921, . 91 4209755704 ,. 9 1 1 796032805 , . 909167983091 

118  DATA  .986595704515, .983939293123, . 90 1 348847046 , .898674465694, . 395966249756 

, .393224301 196, .890448723245, . 387639620403, . 884797098431 ,. 381921264343 

120  DATA  .379012226429, . 376070094 1 95, .373094978418, .870086991109, . 867046245516 

, .363972856122, .860866933638, .S57728610000, .854557988365, .351355193105 

130  DATA  .848120344303, . 344853565250, . 34 1 554977437, . 833224705555, . 834862874986 

, .331469612303, .323045045253, . 324589302785 ,. 82 1 1 0251 4991 , .317584813152 

140  DAI  A  .314036329706, .310457193253, .306347553544, .803207531481, .799537269103 

, .795336904609, .792106577303. .738346427627, .784556597156, .730737228572 

150  DATA  .776838465673, . 773010453363, .769103337646  . 765 1 67265622, . 76 1 2023354*4 

, . 757208846506, . 753136799044, .749136394523, .74505778544;. .740951125355 

160  DATA  . 736816563377, . 732654271672, . 723464390448 , . 724247082951 , . 72000250796’ 

, . 715730325234, . 71 1432195745, . 7071O6731 187, . 702754744457, . 698376249409 

170  DATA  .693971460390, .689540544737, . 685083667773, . 680600997795 ,. 676092703575 

, .671553954847, .666999922304, . 662415777590, . 657306693297, . 653 1 72942954 

130  DATA  .643514401022, . 843331543390, . 6 39 1 24444864 , . 6343932S4 l 64 , . 6296382389 1 5 

, .624859488142, . 6200572 1 1 76 3 ,  6 1 523 1 59053 1 , . 6 1 0392306276 ,. 6055 1 1 04 1404 

190  DATA  .600616479334, .595699304492, . 5907597O1S59, . 585797957456, . 53031 3958096 

, . 575308191413, .570730745-337, .565731810794, .560661576187, .555570233020 

200  DATA  .550457972937, .5453249884 72. .540171472730, . 5349976 1 9887 , . 529303624696 

, .524589632678, . 51 9 3559  90S  66, . 5 1 4 102744 193, . 508930142543, . 503538383726 

210  DATA  .493227666973, .492398192230, . 487550 1 60 1 48, . 4921 3 3772079, . 47679923006 3 

, . 471 396736326, . 465^764957*8, . 4605 387 l 0958 , . 455083587126. .449611 329655 

220  DATA  .444122144570, .433616233539, .433093813853, . 427555093430, . 422000270300 

, . 416429560098, .410343171058, . 40524 l 3 1 4005 , . 399624 l 9984* , . 393992040061 


, . 359395036535, . 354 1 6 3525420, . 3484 18680249, . 342660717312, . 336839853392 

240  DATA  . 331 106305760, . 325310292162, . 319502030316, .31 3631740399, . 307849640042 


302805949319, .296150838244, .290284677254, . 2844075372 1 1 


*8519689335 


250  DATA  .272621 355450, .266712757475, . 2607941 17915, . 254865659605 , . 243927605746 
, .242980179903, . 2 3702 3605994 , . 2 3 1 058 1 0323 1 , . 2250839 t 1 360, . 2 1 9 10 1 240 l 57 
260  DATA  .21 3110319916, .2071 1 1 376132, .201 104634842, . 1 9509032201 6, . 189063664150 
, . 133039337955, . 1 770042204 12,. 170961333760. . 1 649 1 3 1 20490, . 158858143334 
270  DATA  . 152797135253, . 1 46730474455, . 1 40658239333, . 1 345807OSS07, . 1284931 10794 
, .  122410675199, . 1 16313630912. . 1 10222207294, . 104 1 2 1 633872, . 9801 7140 3296E- 1 
230  DATA  . 91 908956497 IE- .$57973 123444E-1 , . 796824 3797 1 4£- 1 , . 735645635997E- 1 , . 
6744  39J95637E-1 , . 6 1 3207 36 3022E -  1 , . 55 1 95244 3497E- 1 , . 49067674 32 74E- 1 
290  DATA  . 429 382569 349E - l , . 3630722294 i 4£- 1 , . 3O674S03 1 766E- t , . 2454 1 2285229E -  1 , . 


134067299058E- 1 , . 1 227153S2S57E-1 , . 6 1 35334649 15E-2,  0 

300  READ  C 

310  K-1024  M 

320  FOR  J"0  TO  H- 4 

330  C<J'*CO«J) 

340  NEXT  J 

350  HI -A  4 

36C  N2«Ml*l 
370  N3*H2* 1 
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